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We present our first results of numerical simulations of two skyrmion dynamics using an u- 
meson stabilized effective Lagrangian. We consider skyrmion-skyrmion scattering with a fixed initial 
velocity of f3 — 0.5, for various impact parameters and groomings. The physical picture that emerges 
is surprisingly rich, while consistent with previous results and general conservation laws. We find 
meson radiation, skyrmion scattering out of the scattering plane, orbiting and capture to bound 
states. 
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. In the large Nc or classical limit of QCD, nucleons may be identified as classical solitons of the scalar-isovector 
SU(Nf) pion field. The simplest theory that manifests these solitons is the non-linear sigma model. However the 
solitons of this theory are not stable against collapse. The first attempt to provide a stable theory was by Skyrme 
' (long before QCD) |l]] who introduced a fourth order term. This term does indeed lead to stabilized solitons that are 
called skyrmions, and there is a vast body of work on their properties and on how to quantize them Unfortunately 
the fourth order term introduces numerical instabilities that make complex dynamical calculation nearly impossible 
It is possible to stabilize the non-linear sigma model without the fourth order, or Skyrme term as it is called, 
■ by coupling the baryon current to the w meson field ||. This not only provides stability, but does so in the context 
J> ' of more reasonable physics. It is also possible to stabilize the solitons by introducing a p meson field, with a gauge 
, coupling, or with both the uj and p, but the p adds a great deal to the numerical complications. Hence in this first, 
OA ■ exploratory work, we stabilize with the uj only. We continue to refer to the solitons as skyrmions. 

The Skyrme approach, either in its fourth order or w stabilized form, has much to recommend it as a model of low 
energy strong interaction physics. This low energy or long wave length domain is notoriously difficult to describe in the 
context of standard QCD. The Skyrme approach includes chiral symmetry, baryons (and therefore also antibaryons), 
pions, the one pion exchange potential, and, with quantization, nucleons and deltas. The idea of having nucleons 
arise naturally from an effective theory of pions and vector mesons is especially attractive as a path from QCD to 
such effective Lagrangians is better understood. 

It has been shown that the Skyrme model can give a good account, with very few parameters, of the low energy 
Oh, properties of nucleons In the baryon number two system, the Skyrme approach can describe the principal features 
of the nucleon-nucleon static potential [Q. Most of these problems have also been successfully studied in more 
traditional nuclear physics forms. 

One problem that has not yielded easily to traditional physics approaches and is out of the reach of standard lattice 
QCD is low energy annihilation. This is a problem ideally suited to the Skyrme approach. We have already shown 
that a general picture of post-annihilation dynamics including branching ratios can be obtained from the Skyrme 
approach [^J. The nucleon-antinucleon potential can be as well H. What remains is a full modeling of annihilation 
from start to finish. The only attempts we know of to do that in the Skyrme model had numerical problems associated 
with the Skyrme term |?j |j| . We propose to study skyrmion dynamics using to stabilization, thus avoiding the usual 
numerical problems. As a prelude to studying annihilation, we have studied scattering in the baryon number two 
system. It is that work we report here. 

In the Skyrme model the three components of the pion field, 7ri,7r2,7r3, can be aligned with the spatial directions, 
x,y, z, providing a correspondence between the two spaces. This simple alignment is called a hedgehog. A rotation 
of the pion field with respect to space is called a grooming. The energy of a single free skyrmion is independent of 
grooming, but the interaction between two skyrmions depends critically, as we shall see, on their relative grooming. 
In this paper we present the results of calculations of skyrmion-skyrmion scattering (with lu meson stabilization) at 
medium energy for a variety of impact parameters and groomings. We find rich structure. Some of the channels have 
simple scattering, but some display radiation, scattering out of the scattering plane, orbiting and ultimately capture 
to a bound B = 2 state. These calculations are numerically complex and require considerable computational resources 
and computing time, but they are numerically stable. We know of no other calculations of B — 2 skyrmion scattering 
that show the phenomena we find. Furthermore our success here bodes well for extending the method to annihilation. 

In the next section we briefly review the model we use, presenting both the Lagrangian for the to stabilized non- 
linear sigma model, and our equations of motion. In Section 3 we discuss our numerical strategies and methods. 
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Section 4 presents our results, mostly in graphical form, and Section 5 deals with conclusions and outlook. The reader 
interested only in results can go directly to Sections 4 and 5. 



II. FORMALISM 
A. Model 

Our starting point is the non-linear a model Lagrangian, 

C a = ^tr(^Wt) + ^tr(W-l) , (1) 

where the SU{2) field U is parameterized by the three real pion fields[| {irk}k=i,z = tt, 

U = exp (i (r • ff)) = cos -k + i sin tt (t ■ it) . (2) 
Here, {rk}k=i,3 = t are the Pauli matrices (in flavor space). We identify the baryon current with 

B" = -^e^tTl(tf dJA) {tfdJA) (ifidfiU)] . (3) 

The full Lagrangian also contains the free oj and the interaction part, 

£ = £ CT + ^uvB" - -d^v - d u w») + -m 2 ^ . (4) 

We take / w = 62 MeV following || and g = -^j- fr om VMD where m = 770 MeV is the vector meson mass. This 
gives a skyrmion mass of 791 MeV. 

B. Choice of Dynamical Variables 

The traditional approach to numerical simulations involving skyrmions uses the cartesian decomposition of the 
unitary matrix U, 

U = $I 2 + * • t ; $ = cos7r, * = 7rsin7r (5) 

The quadruplet of real numbers (<&, \&) is constrained by the unitarity condition UW = l2^ ( E >2 + ^- ^ = l, 
also known as the chiral condition. The Lagrangian is usually written in terms of ($, ty) and the chiral condition is 
imposed using a Lagrange multiplier field. The four coordinates of U are similar to the cartesian coordinates of a 
point in R 4 confined to the surface of a hypersphere of unit radius. 

An attractive idea is to use an approach which ensures naturally the unitarity of U at all times. Such a method is 
used successfully in lattice QCD in the context of the hybrid Monte-Carlo algorithm ||. There, the dynamical variables 
are U itself and the Hermitian matrix 14U^ . In the present work, we employ the parameterization of U that follows 
from (||). We will use tt and 7r and their time derivatives 7T,7r as dynamical variables. This implements only in part 
the principle mentioned above, since the unit vector 7r is subject to conditions similar to the four-dimensional vector 
($, ^>). However, one has better geometrical intuition for the former. The connection to (<£>, ^) is straightforward via 
The Lagrangian in terms of tt, tt is 

£° = l^(d^7r + S m 2 7Td^-d^7r)+m 2 J 2 (cosn-l) 
dm = p,e^ af3 sin 2 ir^ l d„irlT}-(d a nxd 7})} 
C = - \d^ v {d»u v - d'O + im^w" . (6) 



To avoid confusion, throughout this paper we reserve the plain letter 7r for the pion field, and use the symbol 7r = 3.1415 . 
for the mathematical constant. 
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C. Equations of Motion 



We wish to obtain the Euler-Lagrange equations for tt,tt. The problem is that tt is a unit vector, tt ■ tt = 1, so 
its components cannot be treated as true coordinates. One may still write down the Euler-Lagrange equations, by 
considering the action 



S = 



d 3 xdtC (ir(x, t), 7r(x, t), oj a (x, t); d fi n(x, t), d fl Tf(x, t), d IM u} a (x, i)) . 



(7) 



The Lagrange equation for tt is obtained by requiring that S be stationary with respect to a small variation 6ir(x, t) 
and the corresponding variations 8(8^) — d^(5n). The variation of S which has to vanish for any Stt, is 



SS 



d xdt 



—Sir — 
8tt d{d )1 



^y<J(0 M 7r)J = j d 3 xdtl^ 



dC 

to- 9 " 



dC 



Stt 



(8) 



Therefore the quantity in curly brackets in the last equality of (g) has to be identically zero, since Sir is completely 
arbitrary. Of course this leads exactly to the usual form of the Euler-Lagrange equations. We may repeat formally 
the same steps for tt, leading to 



SS = 



<! ^ 8 
8tt 



8C 



8(d^fr) 



Stt 



(9) 



Here, partial differentiation with respect to the vector tt or its derivatives gives the vector obtained by differentiating 
with respect to the corresponding components of tt or its derivatives. If we want SS to vanish for any Stt, we do not 
need to require that the quantity in curly brackets vanish identically. This is because Stt is not completely arbitrary. 
Both tt and tt' = tt + Stt have to be unit vectors. The variation of tt ■ tt must vanish, 

= S(tt-tt) = 2 Stt -ft . (10) 

In other words, Stt(x, t) _L tt{x, t) for any (x, t). Therefore, the necessary condition for the stationarity of S is simply 

8C 



d{d^fr) 

The above statement is easily turned into an equation by subtracting the projection of the left hand side onto tt, 



(11) 
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(12) 



The Euler-Lagrange equations we obtain finally are 

1 



d v d v TT 



■ sin 2tt8 u tt ■ 8^tt — mz sin tt 



3fl 



O^TT = Tt (TT • d^Tt) - ^^daTTdH 



3g 



8 V 8 V ^ 



sm7T 



sin 2 ^ . ^ * x 

^d^d a TT(d p TTXTT) 



l<L e ^ a P S m 2 TTd v TT [TT ■ (8 a TT x dgfr) 
OTT 



(13) 



III. COMPUTATION 

Our calculation is based on the equations of motion ([l3|), using (tt, tt, lu^; tt, tt, w m ) as variables. This choice leads to 
two problems which have to be tackled by our discretization scheme. First, we have a coordinate singularity at tt = 
and tt — tt. Here, tt is not defined. The situation is similar to that of angles in polar coordinates when the radius 
vanishes. While the equations of motion are correct for any non-zero tt, in the vicinity of the coordinate singularity 
small changes of tt are translated into very large variations of tt. A scheme based on the equations of motion ( |13] ) 
breaks down around the 'poles' of the hypersphere due to the large discretization errors involved. 
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One way out is to rotate the coordinate system in SU (2) so as to avoid the problem regions. We introduce a new 
field V obtained by rotation with a fixed Uq S SU(2), 

U = U V ; V = exp(ir-a) ; U = expUr-Qj . (14) 

Substituting (|l4| ) into our Lagrangian (Q)(Q), we find that only the pion mass term is modified since everywhere else 
IA is combined with Vft or its derivatives so Uq drops out. The equations of motion can then be derived in an identical 
fashion to (|l3|). We cite them here for the sake of completeness, 

d v d v a = - sin 27T0 M ff • d^a - ml (cos 9 sin a + (6 ■ a) sin 9 cos a) + sin 2 ad^ v [a ■ (d a a X dpa)} 

d u d»a = a (a- d^a) - ^^d.ad^a ^L ! e! a ' af> duU> v d a a (d a x a) - ml ( 9 - a(Q ■ a)) sin 9 sin a 

* sin % J 2 \ / 

= - m 2 ec ^ - -^e^ a/3 sin 2 a9„a [a ■ (9 Q a x ^a)] . (15) 

In our calculation, we first rotate the field to be updated, together with the surrounding fields, so that the correspond- 
ing (cr, a) is comfortably away from the coordinate singularities. Then we apply the discretized equations of motion 
(see below) derived from (|l^) and finally we rotate back the updated fields. We stress that switching to rotated fields 
and ( |l5| ) amounts to a mere changing of the coordinate system. The content of the equations is identical irrespective 
of the coordinate system. The difference is in the discretization which in the vicinity of the 'poles' leads to large errors 
which are avoided in the rotated frame. The choice of 9 is largely arbitrary. For simplicity we choose it so that a is 
always on the 'equator' of the SU(2) hypersphere. 

The second important numerical issue follows from the unit vector nature^Jof ir, if -if = I. This and the corresponding 
constraints on if and 7r, namely ir -ir = 0, ir -ir + it -if — are consistent with the equations of motion but are violated 
by the discretized equations. Therefore they need to be imposed by projecting out the spurious components of ir and 
ir at every step in the update. There are similar conditions for the spatial derivatives, which also have to be taken 
into account. This is of course a rather technical point, but it is a reminder of the fact that our approach does not 
completely conform to the idea of using a minimal set of coordinates. That would be achieved for instance by trading 
7r, which describes a direction in three dimensions, for the two angles that define that direction. That would have 
introduced another set of coordinate singularities to be avoided. The opposite strategy would be to apply the same 
considerations we followed for ir in deriving the Euler-Lagrange equations and imposing the unit length constraint, to 
the four-dimensional unit vector (<£>, ty). That also remains an open possibility, along with using unitary and Hermitian 
matrices as dynamical variables and of course, the traditional path employing Lagrange multipliers. We have not 
found it necessary to employ any of these notions at this time. Besides the fact that it leads to a reasonably stable 
calculation, our choice of variables has the advantage of fairly simple equations of motion. The decomposition of if 
into its length and unit vector follows quite naturally, and leads to great simplification in the much more complicated 
case of including a vector-isovector p meson field. 

The main idea of our numerical scheme is the following. The discretized time evolution for the fields themselves is 
quite straightforward, given the knowledge of their time-derivatives or 'velocities'. The evolution of the velocities is 
the core of our procedure. The velocities are assigned to half-integer timesteps. Thus, one timeslice contains the fields 
at a given time and the velocities half a timestep earlier. The time evolution of the velocities follows from solving the 
equations of motion for the second time derivatives. The latter are written in terms of the retarded velocity (at time 
(t — Ai/2)) and the one at (t + At/2). The rest of the equations of motion contain the local fields and their spatial 
derivatives, (all defined at time t), but they also contain the velocities as it is evident from (|l5|). We approximate 
the velocities at t by the average of their values at (t ± At/ 2). This leads to an implicit equation for the updated 
velocities. We solve the implicit equation iteratively to cubic order in At which exceeds the order of accuracy in which 
it was derived. Finally the new velocities are used to update the fields. 

Use of a flux-conservative form for the equations of motion might have lead to improved accuracy. However, this 
way we can use the known conservation laws, in particular, energy and baryon number conservation, as a check for the 
accuracy of our simulations. The long-term stability of the calculation is largely improved by the feedback mechanism 
built into our implicit scheme. Calculations performed using an explicit scheme involving a second timeslice shifted 



We will continue to use if when referring to the pion field, even though we perform our updates using the rotated pion fields 

(T. 
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with one half timestep give virtually identical results as far as the time evolution of the fields and energy and baryon 
number conservation are concerned. Part of the calculations presented below - those not involving complicated 
physical situations - have been in fact performed with this earlier version of our code. When pushed beyond 25 fm/c 
in time (a few thousand timesteps), these calculations typically develop instabilities which build up fast to destructive 
levels. By contrast, the feedback calculations ran without exception up to 5000 timesteps or more without becoming 
unstable. This is achieved without adding an explicit numerical viscosity term to our equations. 

The calculation based on the algorithm described above was implemented on a three-dimensional grid of points. 
The physical size of the box used for the calculations shown in this paper was 18 x 10 x 10 fermi. We take advantage 
of the four-fold spatial symmetry of the problem jl0],[ll| , therefore our mesh covers only one quadrant of the physical 
box. We used a fixed lattice spacing of 0.1 fm therefore our lattice has 91 x 101 x 51 regular points. In addition to 
these, we have a layer of unequally spaced points on each outside wall of the box simulating an absorbing boundary. 
Thus our full mesh has 101 x 121 x 61 points. An indication of the better intrinsic stability of the model Lagrangian 
employed here is the fact that our calculations are fairly stable out to 50 fm/c in time with timesteps of 0.1 to 
0.4 fm/c, corresponding to a CFL ratio of 0.1 to 0.4. This is larger than in the early works of Verbaarschot et al JlJ] 
(0.05) and of the Caltech group (h]] (0.075 to 0.013), and comparable to the more recent calculations of jyj, which 
employ fourth-order spatial differences, in contrast with our second-order spatial difference scheme. In the absence 
of radiation, the total energy is conserved to within 3% for typically 20 fm. Due to the emission of radiation which 
eventually leaves the box it is harder to assess the degree of energy conservation for those processes which involve 
(quasi)bound states and have to be followed for a longer time. We can estimate the amount of energy being radiated 
by calculating the energy contained in a sphere of given radius (2 fm) around the origin, large enough to contain most 
of the field. The radiated energy first leaves this sphere and then the box. Excluding loss through radiation identified 
in this manner, the total energy is conserved to within 4 % in all the calculations presented below, including the long 
runs (40 ... 60 fm/c) involving bound states. A check for consistency is to reverse the arrow of time at some point in 
the calculation, which should lead back close to the initial state. We performed this check successfully on one case 
with nontrivial dynamics but little radiation. 

We construct the initial state by first numerically solving the field equations in the spherically symmetric ansatz. 
This gives us the radial functions for the spherically symmetric static skyrmion (hedgehog) . We then place a boosted 
hedgehog configuration in the simulation box. The skyrmion is boosted towards its mirror image implemented via 
the boundary conditions at the symmetry walls of the box, which corresponds only to one quadrant of the physical 
region being simulated. One can dial all the relative groomings discussed below as well as the corresponding skyrmion- 
antiskyrmion configurations by simply choosing the appropriate signs in the boundary conditions for the various field 
components |lj] . The problem of the overlap of the tails of the two skyrmions in the initial state is not easy to solve. 
Instead, we chose the initial configuration with a large separation (9 fm) between the two centers, so that the overlap 
becomes truly negligible. 

Our calculations have been performed on clusters of 6 — 12 IBM SP-2 computers. Depending on the number of 
processors and the timestep, one 20 fm/c calculation takes half a day to two days to complete. 

IV. RESULTS 

For skyrmion-skyrmion scattering with finite impact parameter, there are four relative groomings that are distinct. 
The first is no grooming at all. This is the hedgehog-hedgehog channel. The second grooming we consider is a 
relative rotation of 7r about an axis perpendicular to the scattering plane (the plane formed by the incident direction 
and the impact parameter). The third is a relative grooming by 7r about the direction of motion. For zero impact 
parameter this channel is known to be repulsive. |l3f| The fourth grooming consists of a relative rotation by n about 
the direction of the impact parameter. In the limit of zero impact parameter the second and fourth groomings become 
equivalent, they then correspond to a rotation of 7r about an axis normal to the incident direction. In this zero impact 
parameter case, this channel is known to be attractive. |l3| All of our scattering studies are done at a relative velocity 
of P = v/c = 0.5, which corresponds to a center of mass kinetic energy of 230 MeV. In order to keep our study 
relatively small, we have not studied the effect of varying the incident energy. As the main means of presentation we 
chose energy density contour plots. The baryon density plots are very similar to the energy density plots. 

A. Head-on collisions 

Let us begin, for simplicity, with the case of zero impact parameter. For the hedgehog-hedgehog channel (HH) and 
the repulsive channel (rotation by n about the incident direction) symmetry dictates that the scattering can only be 
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exactly forward (the skyrmions passing through on another) or exactly backward (the skyrmions bouncing back off 
each other). We find that in fact the scattering is backward in both the HH and the repulsive channels. 




FIG. 1. Contour plots of the energy density in the xy plane for scattering in the repulsive channel with zero impact parameter. 
The spacing between the contours is 100 MeV/fm 3 . The first contour is at the 5 MeV/fm 3 level. Note that the frames are 
not evenly spaced in time. Here and in all contour plots, the x axis points to the right and the y axis points upwards, unless 
explicitly specified otherwise. The length on both axes is measured in fermi. 

We illustrate this type of scattering in Figure [l] with contour plots of the energy density in the xy plane^J for the 
repulsive channel. Unless otherwise specified we keep the same choice for the energy contour levels for all similar plots 
that will follow, namely, the first contour at 5MeV/fm 3 , and the others equally spaced at 100MeV/fm 3 . The head-on 
scattering in the repulsive channel has nothing surprising. The process is reminiscent of two tennis balls bouncing off 
each other. The skyrmions start compressing as soon as they touch. They slow down, compress and stop, and then 
expand and move off in the direction opposite to the one they came in along. The collision is practically elastic. We 
were unable to detect any energy loss through radiation and the velocities of the topological centers before and after 
collision are practically the same. 

For the hedgehog-hedgehog channel we again find backward scattering with an evolution very similar to that shown 
in Figure |l|, hence we show no figure. 

In the attractive channel (rotation by 7r normal to the incident direction) the skyrmions scatter at 90° along an axis 
perpendicular to the plane formed by the incident direction and the grooming axis. This right angle scattering is well 
known. |l4| It proceeds through the axially symmetric B = 2 configuration (TB). The skyrmions lose their individual 
identity in this process. 



3 Throughout the paper we adopt the following convention. The direction of the initial motion is x, the impact parameter - 
if nonzero - points in the y direction and z is the direction perpendicular to the [initial] plane of motion, xy. For zero impact 
parameter the choice of y and z is of course arbitrary. 
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FIG. 2. Contour plots of the energy density in the xy plane for scattering in the attractive (1) channel with zero impact 
parameter. The spacing between the contours is 100 MeV/fm 3 . The first contour is at the 5 MeV/fm 3 level. Note that the 
frames are not evenly spaced in time. 

In Figure || we show the energy contours for head-on collision in the attractive channel. At the midpoint of this 
process, shown in in the fourth frame of Figure ^, one can clearly see the torus-shaped configuration. It is situated in 
the plane defined by the incoming and outgoing directions and perpendicular to the grooming axis. The bulk of the 
energy density avoids the center of the doughnut as it shifts from the incoming direction to the one perpendicular to 
it. To illustrate this point better we plot the middle frames of Figure in perspective in Figure 
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FIG. 3. Contour plots of the energy density in the xy plane for scattering in the attractive (1) channel with zero impact 
parameter. These are the same plots as frames 2 through 6 of Figure^, at times t — 7.56, 8.08, 8.44, 8.96, 9.48 fm/c. The 
spacing between the contours is 50 MeV/fm 3 . The first contour is at the 100 MeV/fm 3 level. The distances are measured in 
fm. 

In Figure || we plot the coordinates of the center of one skyrmion (defined as the point where the pion field amplitude 
is exactly 7r) versus, time. The only nonzero coordinate is x, initially. Then, after the right-angle scattering, y is the 
only non-zero coordinate. Straight lines indicate uniform motion. That is the case both before and after the collision. 
The slopes of the y line is noticeably smaller than that of the x line before collision, in other words, the outgoing 
velocity is slightly smaller. This is due to a genuine physical process, radiation, rather than to a numerical artifact, 
because there is no decrease in the total energy of the system. In some of the processes we discuss below, there are 



4 We remind the reader that we plot level contours of the total energy density in the median plane of our three-dimensional 
system, as opposed to three-dimensional surfaces of constant energy. 
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stronger examples of slowing down, accompanied by detectable radiation. 
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FIG. 4. Time evolution of the three coordinates of the topological center of one of the skyrmions, for the head-on collision 
in the attractive (1) grooming. 

As we have seen above, the scattering direction is determined by the incident direction and the grooming direction. 
When the grooming direction is normal to the incident direction, a torus is formed in the plane normal to the grooming 
direction. In the presence of a non-zero impact parameter, it matters whether the grooming direction is parallel or 
normal to that impact parameter. In both cases there is a tendency to the formation of a torus normal to the grooming 
direction, but the ultimate evolution of the scattering is different in the two cases. We will refer to the case where 
the grooming is normal to the impact parameter as attractive (1) and attractive (2) when they are parallel. 

Let us now look at scattering in each grooming as a function of impact parameter. In each of the four cases we 
study impact parameters of 0.4 fm, 0.8 fm, 1.6 fm, and 2.8 fm. In all cases we find that the scattering at 2.8 fm is 
"routine" and hence we do not go to larger impact parameter. 

B. Simple scattering: the hedgehog-hedgehog and attractive (1) channels 




-5-4-3-2-101234 

X coordinate 

FIG. 5. Paths of the topological centers of one of the colliding skyrmions in the xy plane in the hedgehog-hedgehog (no 
grooming) case. The x axis points to the right and the y axis points upwards. The length on both axes is measured in fermi. 

We begin with two groomings that lead to simple dynamics. Consider first the hedgehog-hedgehog channel. In 
Figure |{] we plot the trajectories of the topological center of one of the skyrmions for zero impact parameter and 
for each of the four non-zero impact parameters. Here and in the other similar plots, we define the center as the 
point where the norm of the pion field reaches the value of 7? — n. This is a good indicator of the global movement 
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of a skyrmion, especially when the two colliding objects are somewhat separated.^ We cut off the trajectories in 
Figure |g after t = 17 fm/c. We see normal scattering trajectories corresponding to a repulsive interaction. This 
channel is indeed known to be mildly repulsive. For the smallest impact parameter, b — 0.4 fm, the scattering is in 
the backward direction (recall that for 6 = the scattering angle is 180°) and gradually turns forward for increasing 
impact parameter. This is exactly what one might expect, since this channel is the most similar to point-particle 
scattering. The interaction between the skyrmions is central in the HH channel, because it is independent of the 
direction of the relative position vectors of the topological centers. 

The large-angle scattering for small impact parameters probes the interaction of the soft core of the skyrmions. 
We illustrate this in Figure ^ with energy contour plots from the b = 0.4 fm case. Just like in the head-on case, 
the skyrmions compress as they touch. They slow down and then accelerate and proceed in the outgoing direction. 
Internal oscillations of the skyrmions can be observed after the collision, therefore this process is not entirely 
clastic. 
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FIG. 6. Contour plots of the energy density in the xy plane for scattering in the hedgehog-hedgehog channel with impact 
parameter b = 0.4 fm. The spacing between the contours is 100 MeV/fm 3 . The first contour is at the 5 MeV/fm 3 level. The 
length on both axes is measured in fermi. Note that the frames are not evenly spaced in time. 



The other channel that exhibits only simple scattering is the attractive (1) grooming - where the grooming direction 
is perpendicular to the plane of motion. 



5 The topological center does not necessarily coincide with the center of mass. While the latter is insensitive to internal 
oscillations of the skyrmion, the topological center oscillates. 
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FIG. 7. Paths of the topological centers of one of the colliding skyrmions in the attractive case with grooming around the 
direction perpendicular to the plane of initial movement. The x axis points to the right and the y axis points upwards. The 
length on both axes is measured in fermi. 

We have studied the same five impact parameters here. In Figure ^ we show the trajectory of the skyrmion centers 
for each of those impact parameters. For b = and this grooming, the scattering angle is 90°. For large impact 
parameters the scattering is nearly forward. Hence in between, we should see something in between, which is precisely 
what the figure shows. Note that in this channel, in contrast to the HH channel, the trajectories begin their scattering 
by curving toward the scattering center, as one would expect for an attractive channel. The attraction is seen to act 
practically unperturbed in the scattering with b = 2.8 fm. 




FIG. 8. Contour plots of the energy density in the xy plane for scattering in the attractive (1) channel with impact parameter 
b = 0.8 fm. The spacing between the contours is 100 MeV/fm 3 . The first contour is at the 5 MeV/fm 3 level. The length on 
both axes is measured in fermi. Note that the frames are not evenly spaced in time. 

Again we illustrate one scattering process in more detail. In Figure |^ we plot energy contours for the b = 0.8 fm 
scattering. The most important feature is the existence in the third and fourth frames of a configuration reminiscent 
of the doughnut-shaped intermediate state identified in the b = case for this grooming (Figures || and ||). We are 
further away from axial symmetry than in Figure ||, but the relatively large void in the middle is clearly identifiable. 
It is important to point out that even at this mutual distance of 0.8 fm, the attractive interaction is strong enough to 
start forming the B — 2 configuration. In the present case it is ripped apart by momentum. As we shall see that is 
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not always the case in all groomings. Another feature to mention is the presence of internal oscillations. They can be 
clearly seen in the last few frames of Figure |^ deforming the outgoing skyrmions. In the case of b — 1.6 fm trajectory, 
where the attraction is too weak to lead to a B — 2 configuration, the topological centers are attracted toward each 
other but then they bounce back and oscillate transversally. 

C. Orbiting, capture and radiation: the repulsive channel 

For the remaining two groomings, where the grooming direction is in the scattering plane, the dynamics is more 
complex. Let us first consider grooming about the incident direction, referred to as the repulsive channel. For the 
smaller impact parameters b = 0.4 fm and b = 0.8 fm, we find a remarkable type of trajectory. The scattering begins 
as repulsive as we see in Figure [j| but as the skyrmions pass one another, they find themselves groomed by a rotation 
of 7r about an axis normal to the line joining them. This is the most attractive configuration. Hence as they pass one 
another they begin to attract. The two skyrmions now, feeling this attraction, couple and begin to orbit. In order to 
make the motion clearer in Figure ^, we have decreased the spatial scale of our plot and increased the length of time 
shown on the orbiting trajectories, which are cut off at 28.0 fm/c. 

1 i 1 1 1 1 1 1 
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FIG. 9. Paths of the topological centers of one of the colliding skyrmions in the repulsive case. The x axis points to the 
right and the y axis points upwards. The length on both axes is measured in fermi. 

In Figure [n] we show energy density contour plots for b = 0.8 fm to illustrate our point. Notice that the plots 
cover a longer time period than in the previous cases. As the skyrmions orbit, they radiate. This radiation carries 
off energy and the skyrmions couple more strongly. The radiation and coupling last for a long time. The first energy 
level we plot is 5 MeV/fm 3 , comparable to the amplitude of the radiation. The first burst of radiation can be seen 
in the second and third frames of Figure [l^. It appears to consist of the spinning off of a region of high local energy 
density. A later burst is seen in the last frame. These bursts accompany the orbiting movement of the now bound 
skyrmions. 
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FIG. 10. Contour plots of the energy density in the xy plane for scattering in the repulsive channel with impact parameter 
b — 0.8 fm. The spacing between the contours is 100 MeV/fm 3 . The first contour is at the 5 MeV/fm 3 level so that the 
low-amplitude waves corresponding to outgoing radiation can be observed. The length on both axes is measured in fermi. Note 
that the frames are not evenly spaced in time. 



The radiation takes away some of the angular momentum and a significant part of the energy. The total energy 
radiated is greater than the incident kinetic energy, hence the two are now in a bound state. Thus for this choice of 
parameters, (b — 0.8 fm, v — 0.5 c, repulsive grooming) we have the phenomena of orbiting, radiation and capture. 
We assume that this B = 2 system will eventually find its way to axial-symmetric torus shaped bound state at rest. 
To reach that state they will have to radiate more energy and the remaining angular momentum. We have followed 
the orbiting for a time of more than 60 fm/c. We observed continued orbiting with slowly decreasing amplitude, as 
shown in Figure |ll|, and continuous energy loss through bursts radiation. We are not able, numerically, to follow the 
state to the very end. 
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FIG. 11. Time evolution of the three coordinates of the topological center of one of the two skyrmions. They indicate 
orbiting motion in the xy plane with slow damping. 

We now have a bound B = 2 configuration, and therefore should expect the appearance of a torus. The grooming in 
this case is about the x axis. The torus should now appear in the yz plane, and should rotate about the z axis to carry 
the initial angular momentum. Energy contours in the yz plane corresponding to three frames from Figure [l^ are 
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shown in Figure 12. We recognize the familiar pattern of the doughnut-shaped B = 2 bound state. These snapshots 
correspond to situations when the doughnut is aligned with the yz plane. The doughnut is strongly deformed, with 
the two skyrmions clearly distinguishable. This deformation is only slightly alleviated during the 20 fm/c (and a full 
360° rotation) between the first and the last frame in Figure 
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FIG. 12. Contour plots of the energy density in the yz plane for scattering in the repulsive channel with impact parameter 
b = 0.8 fm. The spacing between the contours is 20 MeV/fm 3 . The first contour is at the 100 MeV/fm 3 . The length on both 
axes is measured in fermi. Note that the plots are slightly tilted to indicate the height of the contour lines. 



The appearance of the B — 2 bound state is even clearer in the b = 0.4 fm case. The xy plane energy contours are 
very similar to those in the b = 0.8 fm process. There is again quite some radiation early on, as illustrated in Figure 
[l3| . As can be seen from the trajectories in Figure |[ the attraction is strong in this case (once the skyrmions are past 
the initial repulsion), and the topological centers come very close. 




FIG. 13. Contour plots of the energy density in the xy plane for scattering in the repulsive channel with impact parameter 
b — 0.4 fm. The spacing between the contours is 100 MeV/fm 3 . The first contour is at the 5 MeV/fm 3 level so that the 
low-amplitude waves corresponding to outgoing radiation can be observed. The length on both axes is measured in fermi. Note 
that the frames are not evenly spaced in time. 

The yz plane contour plot in Figure [lj, showing the doughnut just after its formation, exhibits more axial symmetry 
than the corresponding ones in the b — 0.8 fm case. We conclude that we have an example of two skyrmions merging 
into a B = 2 axial symmetric configuration. The angular momentum remaining after the initial radiation burst is 
carried by the rotation of the doughnut around the z axis. The individual skyrmions lose their identity early on and 
the movement of the topological centers is always very close to the symmetry center. Periodically, they move out of 
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the xy plane in the z direction. We interpret this as oscillations of the torus. We assume that eventually the kinetic 
energy and the angular momentum will be radiated away. 
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FIG. 14. Contour plot of the energy density in the yz plane corresponding to the seventh frame in Figure The spacing 
between the contours here is 20 MeV/fm 3 . The first contour is at the 100 MeV/fm 3 . The length on both axes is measured in 
fermi. Notice that the plot is slightly tilted to emphasize the height of the contours. 

In the b = 1.6 fm case the skyrmions first feel some repulsion, which turns to attraction as they pass. This results 
in small persistent transverse oscillations. For b = 2.8 fm there is only a very small attractive interaction. 



D. Scattering out of the plane of motion: the attractive (2) channel 

The scattering in the last remaining channel, with grooming around the direction of the impact parameter, shows 
the most remarkable behavior of all. The corresponding trajectories are shown in Figure [l5| Recall that for zero 
impact parameter, this grooming leads to right angle scattering in a direction normal to the plane formed by the 
incident direction and the grooming axis. That direction is now normal to the scattering plane, which is the one that 
contains both the impact parameter and the incident momenta. Usually one says that there cannot be scattering, for 
finite impact parameter, out of the scattering plane by angular momentum conservation. However, we have meson 
radiation (mostly pion but also some u>) that can carry off angular momentum, albeit inefficiently. 
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FIG. 15. Paths of the topological centers (in the xy plane) of one of the colliding skyrmions in the attractive case with 
grooming around the direction of the impact parameter. For the b = 0.8 fm case, the return of the trajectories to the xy plane 
is not shown. The x axis points to the right and the y axis points upwards. The length on both axes is measured in fermi. 

For the impact parameter of 0.4 fm this is just what happens. This remarkable trajectory can be better understood 
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by studying the cartesian components of the topological center's position as a function of time. We take the x direction 
as the incident one, the impact parameter in the y direction and the normal to the scattering plane in the z direction. 
The behavior of a;, y, and z as functions of time for the topological center of each of the skyrmions in the case of 
b = 0.4 fm is shown in Figure The skyrmions meet, interact, and then go off normal to the scattering plane in 
the z direction. Therefore they disappear from the plot of trajectories in the xy plane, Figure They disappear at 
the point x = y = 0, with their final state trajectory along the z axis having no impact parameter, as it must not, 
by angular momentum conservation. Note that the slope of the z trajectory in Figure 16 is less than that of the x 
trajectory, reflecting the energy carried off by the radiation. 
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FIG. 16. Time evolution of the three coordinates of the topological center of one of the two skyrmions in the attractive (2) 
channel for b = 0.4 fm. 




FIG. 17. Contour plots of the energy density in the xy plane (top row) and for the same times in the yz plane (bottom row) 
for the attractive (2) channel with impact parameter b — 0.4 fm. The spacing between the contours is 100 MeV/fm 3 . The first 
contour is at the 5 MeV/fm 3 level. The length on both axes is measured in fermi. Note that the frames are not evenly spaced 
in time. 



The energy contour plots shown in Figure |l7] show the skyrmions coming together and attempting to form a 
doughnut about a skewed axis in the xy plane and then flying apart in the z direction. Figure shows the last two 
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configurations at higher resolution in energy. This reveals the radiation that carries off the angular momentum in the 
xy plane. 
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FIG. 18. Contour plots of the energy density in the xy plane and for the same times in the yz plane for the last two 
moments shown in Figure ^ (attractive (2) channel, b — 0.4 fm). The spacing between the contours is 1 MeV/fm 3 , from the 
1 MeV/fm 3 level to the 20 MeV/fm 3 level. The higher energy surfaces are omitted. The purpose of this fine energy scale is to 
show the radiation. 

For the smallest non-zero impact parameter we investigated, b = 0.4 fm, there is meson field radiation left behind 
in the scattering plane that carries off the initial angular momentum. When we go to an impact parameter of 0.8 fm, 
things change. The cartesian coordinates for the b = 0.8 fm case are shown in Figure |l9|. As the skyrmions meet 
there is some curvature and then an attempt at uniform motion along z at x = y = 0. Now there is too much angular 
momentum for the field to carry away. The skyrmions try to go off normal to the scattering plane, but they have only 
a brief excursion in that direction while the meson field is radiating. The skyrmions then return to the scattering 
plane, but by now the field has taken off so much energy that they are bound and they begin to orbit in the xy plane, 
alternating with excursions in the z direction, of slowly decreasing amplitude. 
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FIG. 19. Time evolution of the three coordinates of the topological center of one of the two skyrmions. 

In Figure [l^ we have truncated the trajectory of the b = 0.8 fm case at the point where the skyrmions first leave 
the xy plane. Presumably their final state would once again be the static torus. Between 0.4 fm and 0.8 fm there 
must be a critical impact parameter dividing the cases of skyrmions that escape normal to the scattering plane and 
those that are trapped in orbit in that plane. We are investigating this critical impact parameter and the nature of 
the solutions in its vicinity. 
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FIG. 20. Contour plots of the energy density in the xy plane (top row) and in the yz plane (bottom row) for scattering in 
the attractive (2) channel with impact parameter b = 0.8 fm. The spacing between the contours is 100 MeV/fm 3 . The first 
contour is at the 5 MeV/fm 3 level. The length on both axes is measured in fermi. Note that the frames are not evenly spaced 
in time. 




FIG. 21. Contour plots of the energy density in the xy plane and for the same times in the yz plane for two moments of 
the process shown in Figure [2o| (attractive (2) channel, b = 0.8 fm). The spacing between the contours is 1 MeV/fm 3 , from the 
1 MeV/fm 3 level to the 20 MeV/fm 3 level. 

Energy contour plots for the b = 0.8 fm impact parameter are shown in Figure The first few frames are very 
similar to those in Figure [l^, showing the skyrmions approaching, forming a doughnut practically in the yz plane, 
and attempting to escape in the z direction. As we see in Figure [21], there is considerable radiation at this time as the 
field tries to carry off the angular momentum. The field cannot, and the skyrmions come back first along z and then 
through the doughnut into the xy plane, avoiding the center. Note that both the z turn around and the subsequent 
x turn around comes at no more than 3.0 fm separation, as it must. At larger separation, the skyrmions would be 
out of effective force range and not able to come back. When the skyrmions return to the xy plane, they have lost so 
much energy that they are bound. They begin to orbit, and then make another unsuccessful attempt to escape along 
z. As can be seen from the time evolution plot, this sequence of orbiting of the topological centers in the xy plane 
alternating with excursions out of the plane in the z direction continues for a long time. We interpret this as follows. 
The two skyrmions have essentially merged into a B — 2 configuration. The residual angular momentum forces the 
torus to rotate around the z axis. The remainder of the momentum with which the two skyrmions came into this 
configuration (more precisely, its component pointing to the center), while not enough to push the skyrmions out in 
the perpendicular direction, results in oscillations which deform the doughnut as illustrated in Figure |22j. This central 
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component of the momentum allows the skyrmions to escape in the b = 0.4 fm and the central case. 
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FIG. 22. Two snapshots from the late evolution of the two-skyrmion system in the attractive (2) case, for b = 0.8 fm. The 
xy projections shows that the configuration is close to being aligned with the y axis. The corresponding yz projections exhibit 
the "doughnut" structure. The fact that the deformation is once in the z direction, once in the y direction indicates that 
besides spinning, the "doughnut" also oscillates. The energy level contours here are at 20 MeV/fm 3 , starting at 100 MeV/fm 3 
and the view is slightly tilted. 

We are not able to follow this cycle to completion, but believe that eventually the system will radiate enough energy 
and angular momentum to settle into the static, axial-symmetric B = 2 torus, similarly to the small impact parameters 
for repulsive grooming discussed in the previous section. For impact parameters of 1.6 fm and 2.8 fm in this channel, 
the scattering is unremarkable by comparison, as is seen in Figure [15[ For the largest impact parameter we notice the 
medium-range repulsive interaction, just the opposite of the corresponding case in the repulsive grooming discussed 
above. The two groomings are interchanged at the point where the centers cross the x = plane. 

V. CONCLUSION AND OUTLOOK 

The phenomenology of two-skyrmion scattering for various groomings as it emerges from our qualitative study, 
may be summarized as follows. In the absence of grooming , the scattering process is almost elastic at the energy we 
considered. The kinematics is similar to that of point-particles interacting via a central repulsive interaction. This 
is to be expected, since the interaction between two hedgehog configurations is central. We know [Q that at higher 
velocities one expects to excite vibrational modes of the individual skyrmions. We do see modest indications of that. 

The head-on collision in the repulsive channel is quasi-elastic, similar to the hedgehog-hedgehog processes. The 
remaining collisions involving grooming of 180° can be divided into two categories, depending of whether or not the 
impact parameter, b, is small enough for the formation of the B = 2 (torus-shaped) bound state. If the impact 
parameter is large, the collisions are quasi-elastic with a weak attractive or repulsive character depending on the 
grooming. If the impact parameter is small enough, the collision proceeds through a (sometimes deformed) doughnut 
configuration in the plane normal to the grooming axis, even if the final state is not bound. 

The cleanest example is the attractive (1) channel, with the grooming axis normal to the plane of motion. For zero 
impact parameter it is well known |lj] that a torus is formed in the xy plane, the skyrmions lose their identity, and fly 
out at right angles with respect to the incident direction and the grooming axis. For a nonzero impact parameter, a 
deformed torus still appears in the xy plane, and the scattering angle decreases continuously from 90° as b increases. 

In the repulsive case (grooming around the direction of motion, x), the torus is formed close to the yz plane. As 
the two incident skyrmions approach the y axis (that of the impact parameter), they find themselves in an attractive 
configuration, since they are now groomed around an axis (x) perpendicular to the one (y) connecting them. They 
form the doughnut initially in the plane perpendicular to the grooming axis. The configuration carries some of the 
initial angular momentum by rotating around the z axis. The skyrmions do not tend to exit in the perpendicular 
direction because they came in with very little momentum in the plane of the torus. 

In the attractive (1) and (2) cases, the initial momenta have a large component pointing to the center of the torus. 
This momentum is channeled into the perpendicular direction, always leading to scattering in the attractive (1) case. 
In the attractive (2) case, with grooming around the impact parameter direction y, the doughnut is formed close 
to the xz plane. The 90° scattering would therefore happen in the z direction, but this is strongly limited by the 
necessity of shedding the angular momentum around z. Radiation provides a mechanism for this, and for a small 
enough impact parameter the skyrmions can escape in the z direction. Otherwise they go into a rotating doughnut 
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configuration. In this case however, there is radiation and significant momentum in the plane of the torus, which 
torus also exhibits quadrupolc oscillations. 

While our investigation is by no means complete, we did identify a significant number of distinct patterns for 
the outcome of skyrmion-skyrmion collisions. One may define a number of critical configurations which separate the 
different outcomes, even for fixed velocity. The picture might get even richer by sweeping a range of incident velocities. 

The calculations of skyrmion-skyrmion scattering reported in the Section above demonstrate two things. First they 
show that the u> meson stabilizes the Skyrme model and makes it numerically tractable out to long times and through 
complex dynamical situations. Second the calculation shows a rich mix of phenomena. These include capture and 
orbiting, radiation, and excursions out of the scattering plane. Although each of these arises naturally in the model, 
they have not been seen or demonstrated before, nor have most of them been suggested before. They therefore add 
to the rich and often surprising mix of results in the Skyrme model. Since this is a model of low energy QCD at large 
Nc, the new results give insight into aspects of QCD in the non-perturbative long wave length or low energy domain. 
This is a region in which our best hope for insight comes from effective theories like the Skyrme model. 

The success of these calculations also gives us confidence that the method can be carried over, with only simple 
changes, to the annihilation problem. In particular the stability of the calculations, thanks both to to meson sta- 
bilization and a number of numerical strategies, suggests that the annihilation calculation will also be stable. For 
annihilation, meson field (pion and omega) radiation in the final state is all there is and the fact that we clearly see 
radiation in the skyrmion-skyrmion case is reassuring. 

Our work here suggests a number of further avenues. We have already discussed annihilation. It would also be 
interesting to explore the landscape of skyrmion-skyrmion scattering as a function of energy as well as grooming 
and impact parameter and also to examine more closely the behavior of the model in the neighborhood of critical 
parameters where the scattering behavior changes abruptly. Also interesting would be to try to extract information 
about nucleon-nucleon scattering from the results for skyrmion-skyrmion scattering. We are investigating these 
questions. 
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